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ABSTRACT 

Recent observations and simulations have suggested that Hn regions around massive stars may 
vary in their size and emitted flux on timescales short enough to be observed. This variability can 
have a number of causes, ranging from environmental causes to variability of the ionizing source 
itself. We explore the latter possibility by considering the pre-main-sequence evolution of massive 
protostars and conducting numerical simulations with ionizing radiation feedback using the FLASH 
AMR hydrodynamics code. We investigate three different models: a simple ZAMS model, a self- 
consistent one-zone model by Offner et al. (2009), and a model fit to the tracks computed by Hosokawa 
& Omukai (2009). The protostellar models show that hypercompact Hn regions around massive, 
isolated protostars collapse or shrink from diameters of 80 or 300 AU, depending on the model choice, 
down to near absence during the swelling of stellar radius that accompanies the protostar's transition 
from a convective to a radiative internal structure. This occurs on timescales as short as ^3000 years. 
Subject headings: Hydrodynamics — ISM: HII regions — Radiative transfer — Stars: formation — 
Stars: pre-main sequence 
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1. INTRODUCTION 

Massive stars dominate their birth environments 
through a set of powerful feedback processes, the full 
implications of which are not completely understood. 
Among their various feedback mechanisms is the for- 
mation of Hn regions — bubbles of hot (10 4 K), ion- 
ized gas that expand into the colder (10 2 K) surround- 
ing medium. These regions are fueled by the prodi- 
gious amount s of UV radiation emitted by massive stars 
(jHoare et al1l2007tlBeuther et alj|2007l) . By heating and 
ionizing the gas around them, massive stars alter their 
birth environments and provide detectable observational 
signatures. 

H II regio ns may contribute to the destruction of molec- 
ular clouds (|Ketoll2002L [20ol [2007t lMatzneril2002l ) . help- 
ing to set cloud lifetimes. They are observable by thei r 
radio continuum emission (iMezger fe Henderson! Il967l) 
or by their recombination lines fe.g. lWood fc Churchwelll 
(1989) use the H76a line). Even Young Stellar Ob- 
jects (YSOs) emit eno ugh UV radiation to begin ionizin g 
the gas around them (|Bik et al.ll2005l; IChurchwellll2002D . 
which might teach us something about the lives of mas- 
sive protostars before they reach the main sequence. 
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More recently obser vations have shown time vari- 
ability in Hn regions dFranco-H crnand ez fc Rodriguez! 
1 2004 iRodriguez et all 120071; iGalvan-Madrid et all 120081; 
iGomez et al.l I2008D . iFranco-H crnandcz & Rodriguez 
( 2004) have suggested that such observed time- variability 
may be due to the changes occurring in the source of the 
ionizing radiation, though it may also be due to increased 
absorption in the rapidly-evolving core of the nebula. 

Three-dimensional collapse simulations of massive star 
formation that included feedback by ionizing radiation 
showed time variability with chan ges in size and flux 
very similar to th e observations (peters et al.l I2010at 
Kla ssen et al.l [2012h . The origin of this time variabil- 
ity is the strong accretion flow around the massive star, 
which is dense enough to shield the ionizing radiation. 
As a consequence, the H II region fluctuates between ex- 
tended and trapped states as long as the star is embed- 
ded in a sufficiently strong accretion flow. In the simu- 
lations, the flickering stops when companion stars form 
in the gravitationally unstable accretion flow around the 
central massive star and intercept material that would 
otherwise be accreted onto the high-mass st ar, a process 
we cal l "fragmentation-induced starvation" (Pet ers et al.l 
I2010bf ). Then the accretion flow becomes weak enough 
that the H n region can isolate the forming high-mass star 
and terminate the accretion process entirely. Magnetic 
fields, which provide additional support against gravi- 
tational collapse and reduce the number of companion 
stars formed, do not seem to directly affect the Hn re- 
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gion variability ([Peters et al.ll2011l) . 

The time variability during the accretion phase can 
lead to morphological changes of the H n region appear- 
ance over timescales as short as ~ 10 years and ap- 
pears able to res olve the lifetime problem for ultracom- 
pact H II regions (|Wood fe Churchwelllll989HPeters etafl 
l2010cD . With the aid of a statistical analysis of the syn- 
thetic radio continuum ob servations of the H H r egions 
formed in the simulations, Galvan-Madrid et al.l f|201 lh 
predicted that 7% of ultracompact and hypercompact 
H II regions should have detectable flux increments, and 
3% should have detectable flux decrements when ob- 
served at two epochs separated by about 10 years. Par- 
ticularly interesting a r e shri nking Hn regions, for which 
iGalvan-Madrid et al.l ([201 ID found that flux decrements 
by at least 10% are twice as likely as flux decrements by 
at least 50% for this time interval. 

In general, however, Hn region variability has two pos- 
sible causes. The first is the chaotic gas motions around 
a forming star, where denser, optically thick gas can in- 
tercept stellar radiation and caus ing the shadowed gas 
to neutralize (jPeters et al.l I2010a[ ) . The second possi- 
ble source of H n region variability is due to the star it- 
self, a s already noted bv lFranc~H crnand ez fe Rodriguez! 
( 2004) . Pre-main-sequence stars undergo changes during 
their early e volution, sometimes swelling and other times 
contracting ([Yorke fe Bodenheimerl [20081 ). This can po- 
tentially occur on even short time scales (< 2000 years). 
Changes in the surface temperature of a star will affect 
its UV flux, which in turn affects the size of the Hn 
region. 

The connection between protostellar variability and 
H II region variability has not yet been explored. 
To do so, simulations must be equipped with good 
protostellar models. We ask: can early variabil- 
ity of an Hn region be traced directly to the pre- 
main-sequence evolution of the massive protostar it- 
self? M odels of prestellar evo l ution have been i nvest i- 
gated bvlPalla fe Stahlen (U99l.|Palla fe Stahlerl (fl992l) 
INakano et all (120001) .IMcKee fc Tan! (120031 ) .IQf fner et all 
( 20091) and IHosokawa fc Omukail (12009) , among others, 
but never investigating the H n region variability that 
could result from a mass ive evolving protost ar. 

In a previous paper (jKlassen et all I2012D . we simu- 
lated the effects of using prestellar models in simulations 
of clustere d massive star form ation with initial condi- 
tions as in IPeters et al.l (|2010aD . Fluctuations in H II re- 
gion size and morphology due to the environmental fac- 
tors were unaffected by the choice of prestellar model, 
as might be expected. The prestellar model we imple- 
mented, however, did show a delayed onset of major heat- 
ing and ionization in the cluster. 

In the present paper, we investigate the effect of pre- 
main-sequence evolution on the ionization feedback of a 
massive protostar with radiation-hydrodynamic simula- 
tions that follow the formation and expansion of an H II 
region around a single, isolated star embedded in a ho- 
mogeneous environment. Our reasoning for resorting to 
a simple model calcula tion instead of runnin g a full col- 
lapse simulation as in [Klassen et al. (|2012D is twofold. 
First, it allows us to separate time variability of the 
H II region due to absorption of ionizing radiation by 
the accretion flow from time variability caused by struc- 
tural changes in the protostar. In a realistic simulation, 



both effects would occur simultaneously, making it dif- 
ficult to assess the ir relative importan ce. Second, the 
prestellar model by IQffner et al.l (I2009D . which w as used 
in our collapse simulations (jKlassen et al~ll2012D . repre- 
se nts the swelling phase of the m ore detailed calculations 
by IHosokawa fc Omukail (12009D only very crudely . Fol- 
lowing the behavior in IHosokawa fc Omukail ([2009) more 
closely allows for a better modeling of the time variabil- 
ity. However, accreti on rates in collapse simulations are 
stron gly fluctuating (jPeters et al.l [2010a; Klasse n et al.l 
I2012T) , as opposed to the cons t ant ac cretion rates consid- 
ered by IHosokawa fc Omukail ((2009) . 

Thus, in order to create a realistic model for the time 
variability caused by pre-main-sequence evolution of the 
high-mass protostar, we are forced to use idealized sim- 
ulations with a constant accretion rate as the ones pre- 
sented in this paper. The work is further idealized be- 
cause we deliberately only explore constant accretion 
rates, and thus neglect any time dependence in the HII 
region due to changing accretion. As a result of this sim- 
plicity, we are able to present a clear comparison of three 
sub-grid models for protostellar effects. 

It is at present not possible to run a detailed pre-main- 
sequence evolution code along with a hydrodynamical 
simulation, so that our approach of investigating the two 
aspects of the problem separately appears to be the only 
solution. 

We also stress that we cannot, in general, obtain the 
time dependence of the H n region radius by analytical 
arguments. The key point is that the radius of the pro- 
tostar may grow or contract depending on the stellar 
evolutionary stage, which in turn also depends on the 
accretion history. Since the ionizing flux grows initially 
during the H II region expansion, both an enhanced ion- 
izing flux and the thermal pressure of the already ionized 
gas work together to extend the H II region radius, so that 
the H II region can neither be modeled by a Stromgren 
sphere nor by a standard D-type ionization front. Simi- 
lar considerations hold for the shrinking phase, in which 
the ionized gas partially recombines. 

Generally, a changing stellar radius affects the surface 
temperature, which in turn controls the UV flux. This 
means that if the effective temperature can change sig- 
nificantly over a short period of time, changes in the 
observed flux or size of H II regions may signal evolution- 
ary changes taking place in the stars themselves, whose 
ionizing radiation is driving the region. 

Our numerical approach is described in Section [5J In 
Section |3] we show our results for the evolution of a single 
massive protostar in a uniform medium, accreting mat- 
ter at a constant rate, comparing the effects of different 
accretion rates and different prestellar models. Our as- 
sessment of these effects is discussed in Section [H with 
some considerations of the observability of Hn region 
variability. Our findings are summarized in [SJ 

2. NUMERICAL METHODS 

We perform num erical simulations u sing the flash hy- 
drodynamics code (jFrvxell et al.ll2000l) in its version 2.5, 
which solves the hydrodynamics equations on an adap- 
tive, Eulerian grid. Radiative tr ansfer is handled using 
a raytracing sche me developed bvlRiikhorst et all (2006) 
and extended by IPeters et al.l (|2010aD . We fix a single 
"star" at the centre of the simulation box to serve as the 
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TABLE 1 

Runtime parameters of the single accreting star simulations 



Run Protostcllar Model Accretion Rate Grid Resolution 



la 


Offner et al. (20091 


10" 


3 M /yr 


20 AU 


lah 


Offner et al. (2009) 


10" 


3 M /yr 


10 AU 


lb 


Offner et al. (2009) 


10" 


4 M /yr 


20 AU 


lc 


Offner et al. (2009) 


10" 


5 M /yr 


20 AU 


2i 


Interpolation to ZAMS Table 


10" 


3 M /yr 


20 AU 


3a 


Interpolation to H&O Radius 


10" 


3 M /yr 


20 AU 


3ah 


Interpolation to H&O Radius 


10" 


3 M /yr 


6.5 AU 


3b 


Interpolation to H&O Radius 


10" 


4 M /yr 


20 AU 


3c 


Interpolation to H&O Radius 


10" 


5 M /yr 


20 AU 



equal to the radius of an equivalent-mass star on the main 
sequence. We found that such an approximation for pro- 
tostars will underestimate the ra dius by about an order of 
magnitude (IKlassen et al.ll2012f ). Consequently, estima- 
tions of the surface temperature and ionizing luminosity 
will be strongly overestimated. There is also no variabil- 
ity in the radius besides a steady increase in size as the 
star accretes matter. No structural transitions in the star 
are accounted for and the radius never contracts. To im- 
plement this model in our simulations, we use tabulated 
values of the mass, luminosity and effective temperature 
for a main suquence star. We interpolate within the ta- 
ble based on the current mass of the star and compute 
the radius from the luminosity and temperature. 

The third model we emplo y is based on the results from 
IHosokawa fe Omukail ()2009D (henceforth H&O). In their 
paper, they investigate the evolution of a protostar un- 
der high accretion rates through detailed stellar structure 
calculatio ns. By contrast, rat her than computing stellar 
structure lOffner et all (|2009fl consider the total energy 
evolution of a poly trope. One can c alibra te a one-zone 
model to the IHosokawa fe Omukail fj2009f ) calculations, 
as they describe in Appendix C. However, the accretion 
rates that we test correspond to the published tracks in 
H&O, so we allow our protostellar radius to follow the 
appropriate track using a fitting function. 

2.2. Initial conditions 

The advantage of the radiative feedback code we em- 
ploy lies in its ability to produce realistic Hn regions. 
In a cluster environment, it is difficult to measure the 
effect of the choice of model on H n region size because 
the environment is highly tumultuou s even without any 
initial turbulence or magnetic fields (Pet ers et alj [2010a; 
IKlassen et all 12011) . The H II regions formed in such 
an environment are highly amorphous and time- variable. 
Here, we want to isolate the effects of the protostellar 
model to gauge its impact on H n region formation. To 
explore this interesting possibility, we select a simulation 
box of side length 0.1 pc filled with uniform, isothermal 
molecular gas at 10K and density p = 1 x 10~ 17 g cm" 3 . 
This is not a collapse calculation and self-gravity of the 
gas is switched off. We do not consider any initial tur- 
bulence, nor magnetic fields, nor do we include radiation 
pressure. 

The radiative fe edback code , which is described in 
greater detail in iPeters et al.l (|2010aD . is a hybrid- 
characteristics raytracing scheme that computes gas 
heating with grey atmosphere assumptions, and gas ion- 
ization, based on the computed flux of ionizing photons 
from stars modeled as blackbodies. 

A single star is fixed at the centre of this simulation 
volume whose mass is set to increase at the artificially 
imposed accretion rate M, which is held fixed through- 
out the simulation and is treated as a parameter. By 
artificially fixing the accretion rate in this way, we can 
directly connect H II region variability with the evolution 
of stellar properties, rather than have it be due to fluctu- 
ations in a ccretion rate or envir onment al changes as was 
the case in IKlassen et all (|2012l ). 

We ran simulations with various prestellar models and 
accretion rates in order to compare their effect on H II re- 
gion formation. The various simulations are summarized 



source of radiation and evolve a set of stellar parame- 
ters based on the models described in Section |2~T1 while 
holding the accretion rate fixed. 

The protostellar model provides the radius and the lu- 
minosity, from which we can obtain the surface temper- 
ature of the star, and hence the fraction of photons with 
energies capable of ionizing the surrounding the gas. This 
is not an analytical calculation, but a time-dependent 
simulation that allows us to compute the flux of ionizing 
photons, which the radiative feedback code then uses to 
heat and ionize the gas within the simulation. This al- 
lows for a self-consistent calculation of the size of the H n 
region and its co-evolution with the protostar. 

We perform a grid of simulations with different ac- 
cretion rates and protostellar models. The three pro- 
tostellar models that we selected were (1) a purely 
zero-age main sequence (ZAMS) estimate based on a 
precomputed tabl e of ZAMS v alues, (2) a one-zone 
m odel based on Offner et al.l ()2009D and also used 
in IKlassen et al.l ( 2012D . and (3) a fit to the results 



fro m detailed stellar structure calculations performed 
bv IHosokawa fe Omukail (|2009l ). In all but the ZAMS 
case, we tested accretion rates of M = 10 -3 , 10~ 4 , 
and 10~ 5 M /yr. For the ZAMS case, we only tested 
M = 10 -3 M /yr. For the other two models we have in- 
cluded simulations of M = 10 -5 M Q /yr for comparison, 
even though rates highe r than 10~ Mp, /yr are nece ssary 
to form a massive star (jHosokawa fc Omukail 12009). Ta- 
ble [1] summarizes our simulations. 

2.1. Protostellar models 

In IKlassen et al.l (|2012T) we described the implementa- 
tion of a protostellar evolu tion module for flash based 
on the lOffner et al.l (|2009f) model implemented in the 
ORION code, which we simply called the "Evolving Pro- 
tostar" model. The module considers protostars of mass 
exceeding 0.1 Mq and evolves their stellar parameters 
sclf-consistently with the simulation. The most impor- 
tant of these parameters is the stellar radius. During 
the pre-main-sequence evolution of a massive protostar, 
there comes a point when the stellar structure transi- 
tions from a convective core to a radiative core. This 
is accompan i ed by a swelling of the stellar radius. The 
lOffner et al.l ((2009) model estimates a swelling factor of 
2. 

We compare our result s to a ZAMS model fo r the pro- 
tostellar radius used by IPeters et al.l (|2010a|) . In this 
model, it is assumed that the radius of the protostar is 
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Offner et al. 



Hosokawa & Omukai 
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Fig. 1. — Results of a suite of simulations involving a single, isolated protostar accreting mass at either 10~ 3 , 10~ 4 , or 10- 5 M Q /yr, and 
evolving in stell ar radius (t op row) according to one of three models: an equal-mass ZAMS equivalent (first column), a one-zone mode l 
based on Offner et al. (2009) (centre column), or following the results from the stellar structure calculations by Hosokawa & Omukai (2009) 
(right column). For each of these models, we also show the surface temperature T" e £r, the ionizing luminosity S*, and the radius of the Hll 
region surrounding the protostar. Dotted lines show the ZAMS values. The shaded vertical bars indicate the span in the evolution of the 
protostar between the point where the H II region has reached its maximal pre-collapse size, and the point post-collapse when the previous 
size has been reinstated. In the Hosokawa & Omukai case, the hatched region indicates the span when the Hll region is collapsing. The 
duration of this collapse is indicated at the bottom of the figure in years. In the Offner et al. case, the collapse is instantaneous, so the 
durations at the figure bottom indicate the timescale for regrowth. 
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in Tab le [T] We repeated a portion of th e lOffner et al.1 
(|2009fl and iHosokawa fe Omukail (12009ft models with 
M = lO _3 M0/yr, this time at a high spatial resolution 
of 10 AU and 6.5 AU minimum grid cell size, respectively. 
This better captured the behaviour during the transition 
in stellar structure and swelling of the stellar radius. The 
high-resolution data was then combined with the lower 
resolution data. 

Each simulation was run for long enough to capture 
all the major transitions. In most cases, we followed the 
evolution of the star until it had almost reached the main 
sequence (see Figure [TJ) . 

3. RESULTS 

The data from these simulation runs are presented in 
Figure [IJ which consists of a series of panels, with the 
three columns showing each of the three prestellar models 
we tested. 

The rows of panels in this figure show the key variables 
that we measured or calculated and which describe the 
evolution of the protostar. The first row shows the stellar 
radius in solar units. We see how ZAMS model shows 
nothing but a monotonic increase in stellar radius with 
increasing mass, while the two prestellar models show a 
radius that sometimes swells and sometimes contracts. 
This swelling and contracting of the stellar radius is the 
key to affecting the size of H II regions and we trace the 
effects of this evolution down through the rows of this 
figure. For comparison, in the first row we show the 
ZAMS radius in each panel using a thick dashed line. 

In the second row, we show calculations of the effective 
surface temperature of the protostar. This is computed 
using the equation 

/ L- \ 1/4 

which depends on the intrinsic luminosity of the proto- 
star £ m t and the stellar radius R. The luminosity in the 
first case is estimated as in lOffner et all (p009ft as 

L int = max(L ms , 4tt(tR 2 T^) , (2) 

that is, the greater of either the main sequence luminos- 
ity or the luminosity of a Hayashi-track star with tem- 
perature Tjj = 3000 K. We use this same estimate of 
luminosity for our H&O comparison. The intrinsic lumi- 
nosity in the ZAMS case is calculated by interpolation 
to tabulated values. 

These panels show that when the radius of a star swells 
during the transition to a radiative core structure, the 
surface temperature drops by several thousand Kelvin. A 
ZAMS model misses this effect entirely. Capturing this 
effect in the surface temperature is crucial in correctly 
determining how the ionizing flux S* , shown in row 3, 
changes during pre-main-sequence evolution. 

The ionizing flux of photons from the protostar is cal- 
culated by integrating a blackbody for frequencies above 
the ionization threshold of 13.6 eV to find the number 
of ionizing photons per second. This number depends on 
the surface temperature of the protostar. The prestellar 
models have drops in their surface temperature, which 
correspond to drops in ionizing flux. In the case of 
M = I0- 3 M Q /yr, the lOffner et al.l ([20091 ) model has a 



flux decrement of just over 3 orders of magnitude (from 
5.3 x I0 42 to 1.7 x 10 39 s" 1 ), while the H&O model shows 
a flux decrement of just over 4 orders of magnitude (from 
1.6 x 10 44 to 1.2 x 10 40 s _1 ). The other main difference is 
that in the former case the drop is instantaneous, whereas 
the H&O model show this taking place over 4000 years. 

The ZAMS model has only a monotonically increasing 
ionizing flux and likely also overestimates the ionizing 
flux for much of the pre-main-sequence lifetime of the 
star. 

Finally, we measured the size of the H n region by eval- 
uating the ionization fraction along a line starting at the 
centre of the simulation box, where the particle is lo- 
cated, out along one of the axes to the edge of the sim- 
ulation box. Sampling the ionization fraction along this 
axis, we can determine the radius at which the gas tran- 
sitions from fully ionized to neutral. We define the radius 
of the H II region as the point where the ionization frac- 
tion X is 0.5, which is found using an inerpolation of the 
underlying sampling points. 

The fourth row of panels in Figure [T] shows the key 
result that ties together all the previous findings. Here 
we show the evolving radius of the H n regi on that forms 
around protostar. The prestellar models of lOffner et al.1 
(2009) and H&O, for accretion rates M > lO _4 M /yr, 
show an H II region that collapses or shrinks when the 
ionizing flux is diminished. The loss of ionizing flux is 
ultimately tied back to the stellar radius, which swells 
during the pre-main-sequence evolution of the star. For 
an accretion rate of M = 10 -5 Afo/yr, there was no col- 
lapse of the Hn region. In the H&O case, the growth of 
the H II region at this accretion rate appears stagnated 
for about 100 kyr before it continues to grow. 

The temporary loss of the H n region is more dramatic 
for the higher accretion rates. The region disappears 
completely, falling to near zero from a radius of 150 AU 
in the Offner model, then quickly rebuilding. In the H&O 
model, the H II radius falls to near zero from 40 AU and 
remains undetectable for about 1700 years. 

We added shaded vertical bars that run through Figure 
[T] and mark the key transitions in the prestellar mod- 
els for accretion rates greater than 10 _4 M Q /yr. The 
lengths of time that each pair corresponds to are shown 
at the very bottom of the figure. In the lOffner et al.l 
(2009) model, the shaded regions mark the timescale 
of Hll region collapse and recovery to its pre-collapse 
size. This destruction and regrowth takes about 1700 
years for M = lO _3 M0/yr, and about 6000 years for 
M = 10- 4 M Q /yr. 

In the H&O runs, the shaded region begins when the 
H II region has reached its maximal pre-collapse size and 
ends at the point of full recovery. The first part of the 
shaded region is hatched and ends where the H II region 
has shrunk to zero. The time taken for the region to 
disappear is indicated at the bottom of the figure, and is 
about 3400 years for M = lO~ 3 M /yr, and about 16500 
years for M = 10 _4 M Q /yr. 

In all our model simulations, as the protostar under- 
goes its final contraction towards the main sequence, the 
H II region grows steadily until we terminate the simula- 
tion. 

In Figure [2] we show a closer view of the collapse 
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Fig. 2. — A blown-up view of the time of collapse of the Hn region (bottom row) around the protostar. Left panels: the one-zone model 
by Offner ct al. (2009); right panels: the Hosokawa & Omukai (2009) data, both at an accretion rate of 10 — 3 Mq/yt. The top row shows 
the concurrent effective surface temperature of the protostar. The small-scale variability of i?jj n is partially a numerical effect caused 
when the H II region radius is measured in the post-processing of the simulation data and partially a result of small density perturbations. 
It is on the order of the grid resolution. 



of an H II r e gion around a massive protostar for the 
lOffner et all (|2009f ) and H&O models. This figure in- 
cludes high-resolution data, with minimum grid sizes of 
about 10AU and 6.5 AU, respectively. Time resolution 
was also higher, at dt = 8 years and dt = 7 years, re- 
spectively. In the Offner et al. model, the Hll region 
radius reaches a maximal value of 161 AU, while in the 
H&O case, it grows to only about 45 AU. We show the 
concurrent effective surface temperatures in the top row, 
which shows the H&O reducing in temperature gradu- 
ally instead of instantaneously. When the flux of ionizing 
photons falls due to decreasing surface temperature, the 
H II region shrinks away to zero over a period of about 
3400 years. Only after the surface temperature begins to 
heat up again does the Hll region begin to reform and 
gr ow, as also seen in F igured] The one-zone model based 
on lOffner et al.l (|2009f) has an instantaneous drop in ef- 
fective surface temperature, and also an instantaneous 
drop in Hn region radius. 

4. DISCUSSION 

Of the models we tested, the IHosokawa fc Omukail 
(2009) model is likely the most realistic, since in this 
one we fit functions to their detailed stellar structure 
calculations. The model also contains no discon tinuous 
jumps in stellar radius, as the lOffner et al.l (|2009D model 
does. The timescale for flickering, then, in the nascent 
H H re gion is probably closer to the Hoso kawa fc Omukail 
(2009) results. 

Given the tumulteous environments in which massive 
stars form, variability in the Hll regions seems reason- 



able. IGalvan-Madrid et alJ (|2008l ) reported a substan- 
tial (~45%) decrease in the radio-continuum flux (6 cm) 
of an unresolved H n region around a young massive 
star over a 5 year period. An a nalysis of radiation- 
hy drodynamic simulation data fr om lPeters et al.l (|2010al) 
by IGalvan-Madrid et al.l (|2011l ) using synthetic radio- 
continuum observations showed that H n regions can be 
highly variable over timescales of 10 to 10 4 years. Hn 
regions initially larger than 1000 AU around a massive 
star lost up to 90% of their size in a short span of time 
(~20 years) during a large accretion event. 

Those results indicated that the evolution of ul- 
tracompact and hypercompact H II regions may be 
more complex that the textbook picture of a spheri- 
cal ionized regions expanding into a quiescent medium. 
H II variability is most likely not due to protostel- 
lar evolution, which was sugge s ted a s a possibility in 
iFranco-Hernandez fc Rodrig uez (2004]). The variability 
we observe in our simulations is of a scale too small to 
explain those effects. 

We confirmed that tumultuous birth environments give 
rise to massi ve stars with variabl e H n regions in our own 
simulations (Klass en et al.l l2012). but also experimented 
with a different protostellar model that better captured 
the evolution of the protostellar radius under conditions 
of radidly varying accretion. In this paper, we wanted to 
revisit the toy model of H II regions to see to which de- 
gree a realistic treatment of pre-main-sequence evolution 
affects the H II region formation and evolution. 

We decoupled the environmental effects of the molec- 
ular gas flow through a turbulent environment during 
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the accretion process, to see whether flickering of the 
Hn region could also be a purely evolutionary feature 
of the accreting star. We saw the Hn region around 
the protostar collapse in a similar way, losing over 90% 
of its size i n a m atter of less than 10 years for the 
Offn er et all ff2009h mode l, or over about 3.5 kyrs for 
the ([Hosokaw a fe Omu kai 2009) model. 

Some earlier discussi ons, based on hi ghly symmetric, 
simplified models (e.g. Walmslevl 119951) suggested that 
high accretion rates might "choke off" Hll region for- 
mation inside "hot cores." These hot cores are compact 
(diameters < 0.1 pc), dense (n»-„ > 10 7 c m ~ 3 ), and at 
tempe ratures > 100K (jKurtz et all I2000D . iChurchwelll 
( 2002) describes these regions as the precursors to ultra- 
compact Hn regions, but also suggests that rapid accre- 
tion may prevent young H n regions from begin detected. 
However, it is now clear tha t massive protostars ac- 
crete asymmetricall y via disks (jZinnecker fc Yorkell2007t 
iBeuther et al.ll2007T) . with some fraction of this infalling 
material being launched into ou tflow cavities created 
by jets (jBaneriee fc Pudritzl I2007D . The jet intensity 
will be proportional to the infall rate and will also be 
accompanied by strong beaming of the radiation field 
aroun d the star in what has been called the "flashlight" 
effect (lYorke fe Bodenheimerlll999t lYorke fe Sonnhalterl 
120021: iKrumholz et all I2005D . Analytic models show 
that an Hn region can be sustained in side a bipolar 
outflow cavity despite heavy accretion (Ta n fe McKed 
120031 ) . Numerical simulations of non-s pheric al accre- 
tion (|Peters et all l2010al iKlassen et all 120121 ) demon- 
strate that Hn regions can expand in spite of strong 
accretion flows. 

Under circumstances such as these, the H n region 
would not be choked off, despite the very high accretion 
rate. To our knowledge, accretion choking out the for- 
mation of an H II region has never been seen in numerical 
simulations, nor is accretion ever spherically symmetric 
onto a massive star under realistic conditions. Therefore, 
it is timely to explore the birth and evolution of hyper- 
compact H II regions for massive protostars at high ac- 
cretion rates, even if our spherically symmetric model is 
highly idealized. We sought to explore the time-evolution 
of an H II region around a source whose flux of ionizing 
photons is varying according the time-evolution of the 
protostellar radius and surface temperature. No analyt- 
ical model for the time- variability of H n regions exists 
and that time- variability exists whether the morphology 
is spherically symmetric or not. 

Many numerical calculation s of protostellar evo- 
lution dPalla fe Stahlerl 119911: lYorke fe Bodenheimerl 
200 81 iHosokawa fc Omukail 120091), even for first stars 
(jTan fc McKed I2004D . have shown that the protostel- 
lar radius undergoes a swelling phase followed by a con- 
tracting pha s e as t he star approaches the main sequence. 
iSmith et al.l J2012) considered the time-evolution of the 
protostellar radius for first stars and showed that under 
some conditions the radius can undergo multiple occa- 
sions of swelling and contraction on account of a variable 
accretion rate. We measured the impact of this evolu- 
tion on the size of the nascent Hn region and showed 
that fluctuation in the size does indeed occur, but are 
small enough that they might be very difficult to observe 
with present radio telescopes. 



4.1. Observational considerations 

Hypercompact H n regions are nominally defined with 
sizes of about 30 mpc (6000 AU) or smaller. Understand- 
ing their pl ace in the sch eme of massive star formation is 
important (|Kurtzl 12005). Hn regions, generally, are as- 
sociated with star cluste r s, and we discussed their vari- 
ability in Klasscn et al.1 (|2012| ). Higher- resolution sur- 
veys will ultimately yield more examples of extremely 
compact H II regions associated with individual, massive 
stars instead of clusters that may tell us more about the 
dynamics and kinematics of the regions that give rise to 
a massive star. 

Tan & McKec (2003) proposed that hypercompact H n 
regions would be confined to the outflow cavities cre- 
ated by rapidly-accreting, massive protostars. The faster 
the acccretion, the higher the outflow rat e. In our self- 
consi stent collapsing cluster simulations (jKlassen et al.l 
12012ft . we saw accretion rates that at times approached 
lO~ 3 M /yr. IHosokawa fc Omukail ()2009D also performed 
their calculations of protostellar evolution with accretion 
rates up to 10 _3 M Q /yr. 

Our idealized model of a single star accreting at a fixed 
rate and feeding back into a homogeneous medium iso- 
lates the massive pro tostar from the cluste r environment 
that was studied in IKlassen et al.l (|2012|) . In this iso- 
lated environment, we could see the variability of the 
Hn region, as caused by the protostellar evolution in- 
stead of the environment. It was for this reason that the 
environment was made as simple as possible. Though 
massive protostars in nature accrete material asymmet- 
rically through a disk, this does not affect the fact that 
H II regions will be variable according to their protostel- 
lar evolution, which involves a rapid swelling of the stel- 
lar radius that drastically lowers the surface temperature 
and fraction of photos above the ionization threshold of 
13.6 eV. This also makes accretion onto the protostar 
easier. 

Under the models we tested, the diameter of H II re- 
gions collapsed at rate of 45 AU/yr (for the Offner et 
al. 2009 model), and at a rate of 0.03 AU/yr (for the 
Hosokawa & Omukai 2009 model), where in both cases 
the protostar was accreting at a rate of 10 _3 M Q /yr. The 
Orion Molecular Cloud is located about 460 pc away. 
In observing the radio continuum at 1.3 cm, the EVLA 
has an angular resolution of about 8.7 mas, which cor- 
responds to a spatial resolution of about 4 AU at 460 
pc, and 17.4 AU at 2 kpc. This suggests that Hn region 
flickering aro und massive protos tars could be observable, 
following the lOffner et al.l (|2009l ) model. The Hn regions 
would not appear quite as they do in our simulations, 
given that massive protostars are embedded in dense gas 
and dusk, but they might be seen via outflows as sug- 
gested earlier. 

Additionally, masers serve as useful probes of high- 
mass star-forming environment s and can be as sociated 
with ultracompact Hll regions (jFish et al.l 120051 ). VLBI 
techniques may offer the angular resolution at a few kilo- 
parsecs that ALMA cannot. Multi-epoch observations of 
the hydroxyl masers around W3(OH) with time baseline 
of 8 years showed evidence consistent with an ex panding 
ultracompact Hn region ([Bloemhof et al.lll992f) . Even 
with archival data, the authors were capable of position- 
ing the sources with sub-milliarcsecond resolution, corre- 
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sponding to maser positioning to within better than 2 AU 
at 2.2 kpc. If masers drift is truly associated with Hn 
region variability and could be detected in the earliest 
stages of massive star formation, variability on timescales 
shorter than 10 years could potentially be studied. 

5. CONCLUSIONS 

Our results show that massive protostars in uniform 
media emit ionizing radiation capable of forming a hy- 
percompact H n region. In our simulations using the 
flash astrophysics code with ionizing feedback and dif- 
ferent prestellar models, the pre-main-sequence evolution 
of protostars can cause early variability in the size of 
these H n regions. During pre-main-sequence evolution, 
the stellar radius swells dramatically when the star tran- 
sitions to a radiative internal structure. This is accom- 
panied by a cooling of the surface temperature and a 
drop in the ionizing flux by several orders of magnitude. 
During this transition, young H n regions will shrink or 
collapse until the contracting star has regained its previ- 
ous ionizing flux. 

Our results show that pre-main-sequence evolution 
can affect the growth of H II regions around massive 
protostars on potentially short time scales. The col- 
lapse of the Hil re gion, initially ^300 AU across in the 
Offnc r et al.l (|2009l ) model, happened within a single sim- 
ulation timestep (8 years) when the radius of the star was 
enlarg ed during the evolution. The lHosokawa fc OmukaJ 
(2009) model saw the radius of the Hn region contract 



from an initial diameter of ^90 AU to zero over the 
course of about 3.5 kyr. 

The observability of this variability may be possible 
with the best interferometric radio observations such as 
with VLBI, the EVLA, or ALMA. Nevertheless, model 
simulations of hypercompact Hn regions are an impor- 
tant part of understanding the complete picture of mas- 
sive star formation. Massive protostars are unique in 
that their UV fluxes are large enough to ionize gas even 
before the star has reached the main sequence, and it 
is therefore important that these early states be further 
explored. 
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